# Loading in Packages
library(tidyverse)
library(ggtext)

# loading in User-Defined Functions
source(here::here("my-scripts/R", "eeg_rope_pipeline.R"))
source(here::here("my-scripts/R", "eeg_sim_functions.R"))

Overview of Chapter 1 - Multivarient Approach

Part One - Generate Raw EEG data

In this step, we use 1/F noise generated from a functions derived from Hank Steven’s work (original function is available on GitHub). After define to templates, one containing an effect and one without, we combine these with noise in order to represent a single trial. We then regenerate the noise again for each participant.

Defining Parameters to Create the Noise and Templates

set.seed(314)

true_onset <- 450  ## Effect at 450ms
stim_on <- 300     ## Stimulus is shown at 300ms
freq <- 4         ## Defines resolution of stimulation, frequency of electrode (This is around 250 reading a second)
max_time <- 800    ## End time for Epoch

Xf <- seq(0, max_time, freq) ## start from zero and end at max time, going up in 4ms (freq)
Nf <- length(Xf)   ## Number of timepoints

# Template one only contains zeros
temp1 <- vector(mode = "numeric", length = Nf)  

# Generate ERP Peak
erp <- dnorm(seq(-1.5, 1.5, length.out= 200/freq), 0, 1)
erp <- erp - min(erp)
erp <- erp / max(erp)

# Getting length and storing them for later use
l_erp <- length(erp)
l_pre_stim <- ceiling(true_onset / freq)

# Template 2 contains the ERP peak
temp2 <- c(rep(0, l_pre_stim), erp, rep(0, (Nf - l_erp - l_pre_stim)))

Visualising the ERP Peak

tibble(x = Xf, y = temp2) |> 
  ggplot(aes(x, y)) +
  geom_line()+
  annotate("rect", xmin = true_onset, xmax = (true_onset+200)-freq, ymin = 0, ymax = 1, alpha = .5, fill = "forestgreen")+
  annotate("rect", xmax = true_onset, xmin = 0, ymin = 0, ymax = 1, alpha = .3)+
  annotate("rect", xmax = max_time, xmin = (true_onset+200)-freq, ymin = 0, ymax = 1, alpha = .3)+
  annotate("text", label = "Null Hypothesis", x = 10, y = .95, hjust = 0)+
  annotate("text", label = "Alternative Hypothesis", x = true_onset, y = .95, hjust = 0, colour = "forestgreen")+
  theme_minimal()+
  scale_y_continuous("Time (ms)", expand = c(0,0))

Combining Templates with Noise - Simulating 12 “experiments”, starting with 2 trials and doubling each time (until 4096 trials)

sample_size <-c(10, 25, 50, 100, 150, 200) # number of trials
gsp <- 1 # gamma spectral power 
outvar <- 1 # noise variance

alpha = 0.05 ## Standard Alpha Level for Significance Testing
seed = 1 # set.seed(1)
static_margins = c(0.3, 0.2, 0.1) # static margins for NUll ROPE


## Effect Present in template 2 (temp2)
sim_df_1 <- sim_eeg_rope_pipeline(
                      sample_size = sample_size, # Defines Sample Size/ Number of Trials
                      alpha = alpha, 
                      num_time_points = Nf, 
                      max_time = max_time, 
                      cond1_base = temp1, 
                      cond2_base = temp2, 
                      gamma_spec_power = gsp, 
                      noise_var = outvar,       # Influence Noise output - Zero means no Noise
                      stim_onset = stim_on,     # Time when Stimulus was "shown"
                      seed = seed, 
                      static_margin = static_margins
                      )

## No effect Present (Only using temp1)
sim_df_2 <- sim_eeg_rope_pipeline(
                      sample_size = sample_size, 
                      alpha = alpha, 
                      num_time_points = Nf, 
                      max_time = max_time, 
                      cond1_base = temp1, 
                      cond2_base = temp1, 
                      gamma_spec_power = gsp, 
                      noise_var = outvar,      
                      stim_onset = stim_on, 
                      seed = seed, 
                      static_margin = static_margins
                      )
sim_df_1 |> 
  plot_base()+ 
  labs(subtitle = glue::glue("Lines at Bottom represent significance, 
                            <b span style='color:{line_colours[3]};'>Raw p-values</b span> and 
                            <b span style='color:{line_colours[4]};'>Bonferroni Corrected </b span>"))

Part Two - Define Nulls ROPEs

Experiments with an Effect

We are going to define two ROPE (Regions of Practical Equivalence) using the data to represent the null and alternative hypothesis. The Null ROPE will be defined by the 90% quantile range of the mean pre-stimulus electrical activity for both groups and the Alternative ROPE will be defined by the mean differences of only the significant p-values. It may be better specify the range of time points a priori, however, I will use two different ROPEs to represent the alternative hypthesis, one defined from raw p-values and the other by bonferroni corrected p-values. This is provide a more conservative definition of an effect, as some results pre-stimulus can be significant and those false positive should not influence the alternative ROPE.

Once the ROPE is created for an “experiment”, it will superimpose the data of the next experiment, to see how the new data fits to the previous study. If the Confidence intervals are perfectly contained in a ROPE, then we can consider that as support for rejecting and maintaining the null hypothesis. For example, if the Confidence interval is inside the Null ROPE, even if the difference is statistically significant, this should likely be interpreted as evidence for keeping the null. Likewise, if the alternative ROPE covers the majority of the ERP amplitude, then that would be evidence of a replication of the underlying effect and thus, evidence to reject the null.

sim_df_1 |> 
  plot_nu_rope()
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

# sim_df_1 |> 
#   plot_raw_rope()

As sample size increases the effect becomes more clear and obvious, and the background noise appears to cancel itself out. This in turn increases the size of the alternative ROPEs, both ROPE defined from raw p-values and bonferroni corrected are wided at the larger sample size. However, raw p-values cause the ROPE to cover the entire range of values, including values where an effect is unlikely. Moreover, it appears only the Bonferroni defined ROPE maintains a gap between the two ROPEs, the Null and Alternative, whch would allow easier interpretation. For, if the ROPEs overlap, and the Confidence interval stretches across both, then it is unclear which ROPE it supports, preventing us to make a clear preference of a hypothesis.

It is worth noting that the peak of each of the effects appears to go above the ROPE, as the quantile range which is used to define these ROPEs would need to contain values larger than what is observed. It may not be suitable to use a ROPE to investigate the maximum values, instead it should likely be used to prevent false positives and ensure the observed difference is outside an interval which represents background noise.

Experiments without an Effect

To see how a ROPE can prevent false positives, we ran the EEG simulation but with no underlying effect, therefore any statistically significant results are false positives.

# sim_df_2 |> 
#   plot_bonf_rope()
# sim_df_2|> 
#   plot_raw_rope()

Just like the ROPEs for when there is an underlying effect, the ROPE constructed from raw p-values appears to be less precise and helpful. Here we can see how False positves influence the ROPE whilst no results were significant after bonferroni correction.

Evaluating the False Positives with respect to no ROPE, and the different types of ROPEs

The use of a ROPE should reduce the number of false positives, as if the Confidence interval of a statistically significant point is not completely outside the Null ROPE then no conclusions can be made.

It is probably worth replicating this pipeline (the above process) several times to assess behaviour in the long run - as consistency is difficult to assess from a few trials. These following graphs are influenced by randomness, as the noise is randomly generated, so repeating this pipeline should reduce this unwanted influence.

## Useful to define region where effect is, to know significant results outside that region are caused from noise
effect_time <- c(true_onset - stim_on, true_onset + 200 - freq - stim_on)
seeds = 1:1000
margin = 0.3


## Reading the Pipeline Simulations - 1000 repetitions
sim_rope_1 <- read_csv(here::here("sim_data", "sim_egg_pipeline_eff.csv"))
sim_rope_2 <- read_csv(here::here("sim_data", "sim_egg_pipeline_no_eff.csv"))

## Evaluating the ROPEs of Number of False Positives
eval_rope_1 <- sim_rope_1 |> 
  eval_rope(alpha = 0.05, effect_time_vector = effect_time)

eval_rope_2 <- sim_rope_2 |> 
  eval_rope(alpha = 0.05, effect_time_vector = effect_time)
## Functions to create graphs
source(here::here("my-scripts", "R", "eeg_rope_graphs.R"))
eval_rope_1 |> 
  fp_graphs()
## [[1]]

## 
## [[2]]
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

## 
## [[3]]
## Warning: Removed 5 row(s) containing missing values (geom_path).
## Warning: Removed 5 rows containing missing values (geom_point).

## 
## [[4]]
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 5 row(s) containing missing values (geom_path).
## Warning: Removed 5 rows containing missing values (geom_point).

eval_rope_2 |> 
  fp_graphs()
## [[1]]

## 
## [[2]]
## Warning: Transformation introduced infinite values in continuous y-axis
## Transformation introduced infinite values in continuous y-axis

## 
## [[3]]
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## [[4]]
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

Raw p-values are going to generate false positives, which is expected after doing around 400 t-tests. The ROPE does seem to catch some of these, as if the Confidence interval overlaps with the ROPE, then we don’t conclude it is a significant difference and without a decision. However, bonferroni corrected p-values seems to be prefect as controlling for false positives, which is probably because of the multivariate approach used in the analysis. IF the EEG was analysed with the mean amplitude then bonferroni correction will start to resemble the raw p values, and hopefully the ROPE will still be effect.

eval_rope_1 |> 
  tp_graphs()
## [[1]]
## Warning: Removed 3 rows containing missing values (geom_point).

## 
## [[2]]
## Warning: Removed 3 rows containing missing values (geom_point).

## 
## [[3]]

## 
## [[4]]

# auc_df <- eval_rope_df |>
#   inner_join(tp_df_1, by = c("n_trial", "rope"))

## Not as neat as I expected - maybe we need more values or the data is not properly formatted?
# eval_rope_1|>
#   filter(n_trial != 10) |>
#   ggplot(aes(x = fpr, y = tpr, colour = rope))+
#   geom_line()+
#   theme_minimal()

Next steps: Repeat pipeline to get smooth curves and long-term behaviour of the various methods Repeat pipeline but analyse data using mean amplitude instead of t-tests at every time point Add references to Laken’s and Kruschke’s work, tidy up unprofessional language.

Overview of Chapter 2 - Mean Amplitude Approach

Part One - Generate Raw EEG data

In this step, we use 1/F noise generated from a functions derived from Hank Steven’s work (original function is available on GitHub). After define to templates, one containing an effect and one without, we combine these with noise in order to represent a single trial. We then regenerate the noise again for each participant.

Defining Parameters to Create the Noise and Templates

set.seed(314)

true_onset <- 450  ## Effect at 450ms
stim_on <- 300     ## Stimulus is shown at 300ms
freq <- 4         ## Defines resolution of stimulation, frequency of electrode (This is around 250 reading a second)
max_time <- 800    ## End time for Epoch

Xf <- seq(0, max_time, freq) ## start from zero and end at max time, going up in 4ms (freq)
Nf <- length(Xf)   ## Number of timepoints

# Template one only contains zeros
temp1 <- vector(mode = "numeric", length = Nf)  

# Generate ERP Peak
erp <- dnorm(seq(-1.5, 1.5, length.out= 200/freq), 0, 1)
erp <- erp - min(erp)
erp <- erp / max(erp)

# Getting length and storing them for later use
l_erp <- length(erp)
l_pre_stim <- ceiling(true_onset / freq)

# Template 2 contains the ERP peak
temp2 <- c(rep(0, l_pre_stim), erp, rep(0, (Nf - l_erp - l_pre_stim)))

Visualising the ERP Peak

tibble(x = Xf, y = temp2) |> 
  ggplot(aes(x, y)) +
  geom_line()+
  annotate("rect", xmin = true_onset, xmax = (true_onset+200)-freq, ymin = 0, ymax = 1, alpha = .5, fill = "forestgreen")+
  annotate("rect", xmax = true_onset, xmin = 0, ymin = 0, ymax = 1, alpha = .3)+
  annotate("rect", xmax = max_time, xmin = (true_onset+200)-freq, ymin = 0, ymax = 1, alpha = .3)+
  annotate("text", label = "Null Hypothesis", x = 10, y = .95, hjust = 0)+
  annotate("text", label = "Alternative Hypothesis", x = true_onset, y = .95, hjust = 0, colour = "forestgreen")+
  theme_minimal()+
  scale_y_continuous("Time (ms)", expand = c(0,0))

sample_size <-c(10, 25, 50, 100, 150, 200) # number of trials
gsp <- 1 # gamma spectral power 
outvar <- 1 # noise variance

alpha = 0.05 ## Standard Alpha Level for Significance Testing
seed = 1 # set.seed(1)
static_margins = c(0.3, 0.2, 0.1) # static margins for NUll ROPE